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FILTERING ANALYSIS OF A DIRECT NUMERICAL SIMULATION 
OF THE TURBULENT RAYLEIGH-BENARD PROBLEM 


T. M. Eidson 1 
M. Y. Hussaini 1 
T. A. Zang 

ABSTRACT 

A filtering analysis of a turbulent flow has been developed which provides details of 
the path of the kinetic energy of the flow from its creation via thermal production to its 
dissipation. A low-pass spatial filter is used to split the velocity and the temperature field 
into a filtered component (composed mainly of scales larger than a specific size, nominally 
the filter width) and a fluctuation component (scales smaller than a specific size). Variables 
derived from these fields can fall into one of the above two ranges or be composed of a mixture 
of scales dominated by scales near the specific size. The filter is used to split the kinetic 
energy equation into three equations corresponding to the three scale ranges described above. 

The data from a direct simulation of the Rayleigh-Benard problem for conditions where 
the flow is turbulent is used to calculate the individual terms in the three kinetic energy 
equations. This is done for a range of filter widths. These results axe used to study the 
spatial location and the scale range of the thermal energy production, the cascading of 
kinetic energy, the diffusion of kinetic energy and the energy dissipation. These results 
are used also to evaluate two subgrid models typically used in large-eddy simulations of 
turbulence. Subgrid models attempt to model the energy below the filter width that is 
removed by a low-pass filter. 


1 Research was supported by the National Aeronautics and Space Administration under NASA Contract 
No. NAS1-18605 while the authors were in residence at the Institute for Computer Applications in Science 
and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 
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1 Introduction 


The numerical simulation of turbulent flows has provided considerable insight into flow 
physics over the last 20 years[13, 14]. The Rayleigh-Benard problem is an example of a 
specific flow which is better understood because of numerical experiments [7]. The high cost 
of direct numerical simulations (DNS) and the lack of confidence in extrapolating subgrid 
models in large-eddy simulations (LES) to new flows has prevented this technology from 
being more widely used. 

A DNS of a turbulent Rayleigh-Benard flow was computed [8] to provide a database 
for a more extensive analysis of a turbulent flow. The specifics of the simulation and some 
preliminary data analysis have been presented previously [8]. The focus of the analysis pre- 
sented in this paper is to use low-bandpass filtering to study the flow field generated by a 
direct simulation. The intent is to gain insight into the turbulent flow simulated as well as 
to understand better the filtering and modeling used in large-eddy simulations. 

The LES approach is based on splitting the velocity field into a large-scale (filtered 
velocity) and a small-scale (fluctuation about the filtered velocity) component. The large- 
scale component is calculated from the appropriately modified Navier-Stokes (and thermal 
energy) equations (LES equations) and is resolvable. The modification is basically done by 
adding a subgrid model for the effects of the small-scale motion. A rigorous derivation of 
the LES equations proceeds by spacially filtering each term in the original equations with a 
low bandpass filter. Details of the derivation of the LES equations can be found in [14]. The 
result is a set of differential equations for the filtered variables (usually the velocities and 
the temperature) with any terms involving a fluctuation variable replaced by an appropriate 
subgrid model. 

To gain insight into the LES approach the total dependent variable fields can be computed 
using the DNS approach and the terms involving fluctuation quantities in the LES equations 
calculated to evaluate proposed models. This has been previously done for several flows for 
which the fluctuation terms (or “turbulent stresses”) were directly calculated and compared 
to models by a correlation analysis, for example [1, 4], In the present study a correlation 
analysis similar to those of previous studies is done with the following differences. Since the 
vertical direction is inhomogenous, the correlations are only integrated over the horizontal 
directions. Thus the correlation coefficients are functions of the vertical coordinate. Since the 
correlation analysis does not compare the mean value of the terms of interest, the horizontal 
average of these terms are calculated using both a model and the “exact” turbulent stress. 
These horizontal averages are compared directly to evaluate the model[ll]. Also the fine 
grid data is frequently “smoothed” onto a coarse grid with a grid size approximately equal 
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to but smaller than the filter widths used. This is not done in this study in order to separate 
the effect of the numerical method from the effect of the filter width. 

While insight into the LES technique is one motivation for this study, the filtering analysis 
is independent of LES. A range of filter widths from the grid size (or viscous scales, whichever 
is larger) up to the size of the largest flow scales of the problem can be used to analyze the 
data from a particular flow simulation. The terms in the various kinetic energy equations 
based on the filtered and fluctuation velocities can be calculated as functions of the filter 
width. This yields information similar to that of energy spectra is obtained but with a 
different perspective. 

The filter analysis results presented here are calculated from the output of a DNS of a 
turbulent Rayleigh-Benard flow (natural convection) at a Rayleigh number (Ra) of 3.8 x 10 s 
and a Prandtl number (Pr) of 0.76 . These values were chosen because of the availability of 
experimental and other simulation data at similar values. Details of the flow, basic variable 
defintions, and the numerical simulation can be found in Appendix A and in a previous 
publication [8] where the simulation results were directly compared to other available data, 
both experimental and numerical[5, 2, 9]. 

Briefly, the simulation entailed the computation of the 3-dimensional velocity and the 
temperature fields on a 128x64x64 grid. The ratio of horizontal to vertical length for the 
computational volume was 4 and 2 in the x\ and X 2 directions, respectively. After a steady 
state flow was developed, data was collected over a time period equal to 10/ W C) where 
W c = (NuPrRa) 1 / 3 . This period, which should consist of several large eddy turn-over 
times, was found adequate by Eidson[7]. The average Nusselt number (Nu) calculated from 
the simulation data is 6.6 . 

The analysis reported previously [8] shows that the simulation output is well resolved. 
The global quantities, the Nusselt number and the RMS averages of the velocity and temper- 
ature, are in good agreement with the experimental results. The Nu is the closest of several 
LES and DNS studies to the average of several experimental studies (Nu = 6.0). The vari- 
ation in x 3 (the vertical coordinate) of the various velocity and temperature quantities and 
the terms of the total kinetic energy equation are in good agreement with the experiments. 
Most encouraging is the comparison with variance and skewness measurements of the tem- 
perature and its vertical derivative by Carroll[2]. While the agreement in magnitude is only 
fair, the changes in slope with x 3 are excellently predicted. These results along with some 
consistency checks for the volume averaged terms in the various kinetic energy equations 
which are made in the current analysis provide a high degree of confidence in the simulation 
results. 
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2 Derivation of the Filtered Kinetic Energy Equa- 
tions 

The choice of filter to separate the large-scale and small-scale components is important. The 
fundamental structures or basis functions of turbulence are not sufficiently understood to 
select the best filter to separate large-scale from small-scale flow structures. Sharp cut-off 
filters have been used with some success for LES studies, but they do not allow for flow 
scales near the filter width to be treated separately from the smaller scales. Consistent with 
the current understanding of turbulence any smooth low bandpass filter is satisfactory. A 
Gaussian filter function (shown for 1 dimension) is generally used and thus was selected for 
this study. 

u(x) = I G(x — x)u(x)dx (1) 

J — oo 



G is the filter function and is the filter width. The above filter is applied only to the two 
horizontal directions, xi and x 2 , to avoid the complications of filtering in an inhomogeneous 
direction, x 3 . The fluctuating component (denoted by a prime) of a variable, e.g. u(x,t), is 
defined as follows: 

u'(x,t) = u(x,t) — u(x, t) (3) 

Five filter widths are used to analyze the simulation data — 1/16, 1/8, 1/4, 1/2 and 3/4. 
For comparison, the length of the sides of the control volume for the simulation, L Xi , are: 

L XI = 4 

L X3 = 2 

L Xi = 1 (4) 

The grid dimensions in the horizontal directions are both 1/32 . In the vertical direction a 
non-unform Chebyshev grid is used with the average grid size being 1/64 . 

The derivation of the filtered kinetic energy equations proceeds as follows: 

1. The time-dependent partial differential equations for continuity and momentum (Equa- 
tions 78 and 79 in Appendix A) are filtered term by term. 
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2. Filtered velocity products are expanded as follows: 


UiUj = U x Uj + U^Uj + UiUj + U-u'- 


— + ~Lij + TJij + TZij 

( 5 ) 

Rij = u'u' 

(6) 

Uj -j- U^Uj 

( 7 ) 

Mij = UiUj 

(8) 

Zij = Hij - Mij 

( 9 ) 


3. Filtered derivatives of a variable can be shown to be equal to derivatives of the filtered 
variable. 


~5u du 
dx i dxi 

4. Using 2. and 3. to simplify the result of 1., filtered equations are formed. 

dUi 


( 10 ) 


dxi 


= 0 


( 11 ) 


9itj dMjj 
dt dxj 


ap + pM 


dxi dxj 

j + Rjj + Ljj) 

dxj 


+ PrRaTSia 


( 12 ) 


5. The filtered equations are subtracted from the instantaneous equations to form fluctu- 
ation equations. 


£4 = 0 

OXi 


( 13 ) 



dt 


dP' „ ds< 

dxi dxj 


agjj - Cjj) a (ftj - Rjj) dig 
dxj dxj dxj 


+ PrRaT'Si 3 


(14) 
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6. A kinetic energy equation is derived from a dot product of a velocity times the appro- 
priate momentum equation. The total equation is the total velocity times the instan- 
taneous momentum equation. The filtered equation comes from the filtered velocity 
and the filtered momentum equation and the fluctuation equation uses the appropriate 
fluctuation quantities. The cross kinetic energy equation is formed as the fluctuation 
velocity times the filtered momentum equation added to the opposite (filter velocity 
and fluctuation momentum equation). 


• Total Kinetic Energy Equation 

dE _ fl[uj(ig + P)] p r d{uiSn) 
dt dxj dxj 

+PrRa(u 3 T ) - ^PrSijSij 

= D* + P T -e (15) 

E = ^u iUi (16) 


D* 


dME + P) } dKftj) 

dxj dxj 


(17) 


P T = PrRa(u 3 T ) 


(18) 


6 ~ \ PrS 'i Si j ( 19 ) 

• Filtered Kinetic Energy Equation 

dEp_ _ d[uj(i?F + -P)] p r ^jv^ij) 
dt ~ dxj dxj 

+ PrRa ^ 3 f ) 

vlj 

= D* f + W F + P$ - t F 
= D F + P F + P^-€ F 

E f = ^UiUi 

n . F = d[*i(Er + P)3 

dxj dxj 


( 20 ) 

( 21 ) 

( 22 ) 
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( 23 ) 


rxF _ jy*F _ du,(Z7,j; + Rjj + Ljj ) 

dxj 

W F = wg' + W f R } + W/' 

= - P ff _ p// _ p[f 


Pf = PrRa'{u 3 T ) 


«" = 5^3.,- 


Cross Kinetic Energy Equation 

dEc _ dju'jP + U]P') | p r d(ujSjj + UiS<j) 


dt dxj 

ePa-Cn) ac* 
+Ui d 7 j “<■ 

_L 7 7 ^(^tj ~ ^*j) __ / d-Rjj 

+ ‘ dxj ‘ 5x,- 


dxj 


dL; 


+Ui ~^ ~ u 'n£r + PrRa '( u 3 T + «3 t') 
-PrSiM 

= £>* c + W c + P£ - e c 
= D c + P c + P?-e c 


Ec = u\u{ 


, c = a(u ;p + s,f) aw?,, + s^.) 

5sj r 


D* c = 


D c = + appropriate stress diffusion terms 


(24) 

(25) 

(26) 
(27) 


(28) 

(29) 

(30) 

(31) 


W c = (W? - Wl’) + W’J + (W'‘ - If") + W# 
+(Wtf ~ Wl>) 


(32) 


P c = 


-{Pt‘ - Pc) - PS’ - (Pr - PH) - pH 
~(P r J - p[') 


( 33 ) 
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P$ = PrRa'lu'^f + u 3 T') 


(34) 


e c = PrlijS;, (35) 

• Fluctuation Kinetic Energy Equation 

m aw) a(„i5{,) 

dt dxj r dxj 

.spa -cj 

— + “‘ — aT, — 

+<^f- + PrRa'^T') 

= D* s + W s + P*-e s 


= D s + P s + Pf-e s 

(36) 

E S = \u\u\ 

(37) 

jy - a K p ’> , *.#«%) 

dxj ' d Xj 

(38) 

D s = D* s + appropriate stress dif fusion terms 

(39) 

W s = (Wg - Wg) + (Wg - Wg) - Wg 

(40) 

P s = -{Pg - pg) - (Pg - pg) + pg 

(41) 

P$ = +PrRa'(u' 3 r) 

(42) 

t s = \PrSPS', 

(43) 


The production terms are composed of the thermal production, Pt and the velocity 
production, Pg. The velocity production takes the following form: 

Pt = -5 ■ (44) 
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where is one of the “turbulent stresses” — the Reynolds stress, the Cross stress, 
or the Leonard stress, Lij. Another term which appears in the analysis with a role similar 
to a stress is M{j. The superscript b refers to the filtered value of the stress(6 = /) or to the 
unfiltered or total stress(6 = t ). These terms are called production terms since they have 
the role of extracting energy from the flow at one velocity scale and adding energy to the 
flow (or producing energy) at another scale. The term Sij a is the strain tensor which can be 
formed from filtered velocities (a = /), ~Sij, or from fluctuation (or perturbation) velocities 
(o = p), Slj. The notation for each of the particular production terms that appear in the set 
of filtered kinetic energy equations derived in this paper is given in Appendix B. 

The production terms are only one form that the work of stresses can take in a kinetic 
energy equation. Another form is as follows: 

r)7 b 

Wf = (45) 

This form will simply be called the work term. The work term can be expanded into the 
production term plus another work term, T£ 6 which measures the diffusion of the work of 
the stresses. 


Wf = Tf - Pf 


rpab 


z — 


dxj 


(46) 

(47) 


Tf will be referred to as the turbulent stress diffusion. The volume average of the production 
term is the negative of the volume average of the work term for the same Z and (at). Using 
the straight-forward derivations for the kinetic energy equations described above, the work 
form naturally appears. Equation 46 can be used to modify point wise the kinetic energy 
equations to include the production form. In general, the production form and the work 
form are not equal locally. One exception is the following: 


P p r = Wg (48) 

Each of the kinetic energy equations can be written as “the time derivative of a kinetic 
energy equals diffusion terms plus production (or work) terms minus a dissipation term.” 
If the flow has reached a “turbulent steady state” then the time derivative of the kinetic 
energy averaged over a suitably large spacial region should be zero. The diffusion terms are 
of the form — the sum of the spacial derivative in each of the 3 directions of a group of terms. 
The volume average of the diffusion terms can be shown to equal zero for the boundary 
conditions in this simulation; specifically, periodic boundary conditions in the horizontal 
directions and zero velocities at the walls. Otherwise, these terms are not zero even for a 
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horizontal average. The disappearance of the time derivative and diffusion terms for suitable 
averages partially explains why a “production equals dissipation” assumption is often made 
in turbulent theories [12, 10]. All the terms in these kinetic energy equation are evaluated in 
this study. 

The filtered momentum and continuity equations become the LES equations after the 
appropriate terms involving fluctuation variables have been modeled. Actually the fluctua- 
tion kinetic energy equations (and even equations for R+j and C{j) with appropriate models 
can be incorporated into the LES equations, but this has not been generally done[6]. For 
the present study Equation 15 will be called the total equation since the principal dependent 
variable, E , is a function of total velocities, u t . The filtered kinetic energy equation also can 
be called the resolvable scale equation since the principle variable, Ep, is a function of the 
filtered velocities, which would be directly resolved in the LES technique. The perturbation 
kinetic energy equation can be viewed as a subgrid equation since the principle variable, Es, 
is a function of velocities with scales below the filter width and for LES, the filter width is of 
the same order as the grid size. The cross kinetic energy equation can be called the transfer 
scale equation since it has been suggested that Ec is principally composed of velocity scales 
that transfer energy from the filtered region to the fluctuation region[l]. Equations 20, 28 
and 36 will be referred to as the set of filtered kinetic energy equations. 

The total potential energy equation is included for reference. The potential energy equa- 
tion is derived by multiplying the thermal energy equation by —PrRa(x 3 ). Note that the 
term, PrRa(u 3 T T ), removes energy (has a negative sign) from the potential energy equation 
where it adds energy (has a positive sign) to the total kinetic energy equation. 

fls 3$ ( d 2 T \ 

m +u ’aTr - PrR <^ - PrR ‘ (- 57 ) (49) 

$ is the non-dimensional potential energy. 2 

$ = -PvRaT r x 3 (50) 

2 In dimensional form, $ is defined as follows: 

* = = - g/3T r x 3 

Po 
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3 Data Analysis Details 


When the kinetic energy equations are horizontally averaged, the cross kinetic energy equa- 
tion can be simplified. The simplification is a consequence of the following identity. 

< fg > = J f(x) G(x — x')g(x')dx'^j dx 

= J g(x') {^J G{x — x')f{x)dx^j dx' 

= < f 9 > ) (51) 


where f and g are smooth, continuous functions of x. The < > denote a horizontal average. 
When the above identity is applied to certain production and work terms the following 
relationships can be derived. 


< Ui 


dL 




dxj 


> = < Ui- 


dM 




dM, 


= < Ui 


dxj 

mi 

dxj 


- ^ 




dxj 


> 


- Ui- 




dxj 


,SMa 

< U ‘ dx 5 > 


(52) 


e- fr -^»j) ^ 

< u i— a > 

OXj 


dRii dRii 

dxj dxj 

dRij _ dRij 

,d~R,j 

~ <Ui dx j > 


<fli g,,) >=<u\?h> 


dxj 


dxj 


< ~3.il,, >=-< S'.iM.j > 

< - ~R,i) >=< SliHn > 

< >=< SI, Hit > 


(53) 

(54) 

(55) 

(56) 

(57) 


These relationships can be used to simplify the horizontally averaged cross kinetic energy 
equation to the following form: 


dEc 

< ~df 


>=< D* c + 2 Wp + 2Wr + 2W V J + P% - t c > 


(58) 
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or 


< £!££ >=< D c - 2 Pq } - 2 P# - 2 P$ + P$ - e c > (59) 

ot 

Although these equations are derived for continuous arithmetic, the derivation for discrete 
arithmetic is straight-forward. 

The data analysis proceeded as follows: 

1. The velocity and temperature values at each grid point were saved for 322 time steps 
during a “steady state” portion of the DNS. These time steps are equally spaced over 
a time period of 10/W C . For each filter width, the various production, work, diffusion 
and dissipation terms are calculated at each grid point for a particular time step. All 
derivatives and multiplications of variables are done in the same manner as the original 
simulation[8]. The analysis arithmetic is done using 32-bits only compared to 64- bits 
used in the original simulation. 

2. The terms are then horizontally averaged, denoted by < >, and volume averaged, 
denote by { }. The data was then “long” time averaged over about 80% of the “steady 
state” portion of the simulation — the time period shown in Figure 1. 

In the plots of the horizontal averages versus z 3 , the curves are not always symmetrical 
about 13 = 0, particularly the production/work terms (e.g. Figure 16). Specifically, some of 
the variations with x 3 (or “humps”) appear to be too irregular to be a “general steady state” 
solution that would be found regardless of initial conditions (and possibly small changes in 
boundary conditions such as the aspect ratio of the computational volume). The horizontal 
and volume averages of the velocities should tend to zero since viscosity should dissipate any 
mean horizontal motion. “Non-symmetric or non-zero” initial conditions or a weak organized 
flow caused by an instability may persist for a long period of time[3]. Plots of the horizonal 
velocity averages showed that a small patterned flow did exist during the entire simulation. 
In Figure la, {ui} and {it 3 } are small compared to the RMS levels in Figure lb but show no 
tendency to dissipate. A change in the pattern of {u,} was observed in Figure lb between 
the first and second half of the simulation. Differences in < tt; > versus x 3 between the two 
halves are also observed. The pattern in Figure 2a for 4 selected time steps is observed during 
the entire simulation with the peaks oscillating several times from positive to negative. This 
pattern could have resulted from a large cellular motion that filled the entire vertical region. 
The pattern in Figure 2b is observed from the middle to the end of the simulation. This 
form could result from cellular motion of about half the size of the distance between the 
plates. < u 3 > is several orders of magnitude smaller than < > and < u 2 > due to 
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the effect of the continuity equation and the periodic boundary conditions in the horizontal 
directions. These facts lead one to suspect that the component of the flow causing the non- 
zero magnitude of < U\ > and < u 2 > is related to the irregular “humps’ in the velocity 
production terms. < Ui > and < u 2 > are a maximum of about 20 at each time step. The 
production terms are proportional to velocity cubed and the “humps” are generally of order 
20 3 = 8000. Also the x 3 size of the “humps” is consistent with the variations in x 3 seen in 
Figure 2b as well as similar plots at other time steps (not shown). 

One other source of inaccuracy that is investigated is the effect of aliasing. Although 
there is no apparent aliasing problems in the simulation, several of the operations in this 
analysis — multiplication, taking derivatives, extracting perturbation variables — increase the 
relative size of the higher wavenumber components of the field being calculated. As a check, 
the original velocity fields were expanded by spectral methods onto a grid with twice the 
resolution and a limited set of results were calculated on this grid. Although some differences 
in the results are observed locally (a small 3-dimensional spacal region at one time step) for 
the most sensitive quanties being calculated, no significant difference is found between the 
horizontal averages calculated on the original grid and on the expanded grid. 

4 Volume Averages 

The data calculated for the five filter widths specified previously are discussed below where 
the terms of the kinetic energy equations have been averaged over the computation control 
volume. The filter widths, listed in Table 1, can be compared to the size of the control 
volume for reference, which is l(height) by 4 by 2 in units of non-dimensional length. The 
data is long time averaged. 

In Table 1 the non-zero volume averages for the terms in Equations 20, 28 and 36 along 
with the kinetic energy associated with each equation are shown. Also the values for the 
diffusion terms, which should give zero, are included as a measure of the numerical accuracy 
of the original simulation and current analysis. In Table 2 the fraction of kinetic energy in 
the scales below the filter width is given. 

These volume averages are plotted versus filter width in Figures 3a-d. For the small filter 
widths the filtered kinetic energy (and terms in the associated equation ) dominates with the 
cross kinetic energy (and terms) being larger than the fluctuation equivalent. At the larger 
filter widths, the energy is more evenly divided. In each of the 3 kinetic energy equations 
the thermal production and the dissipation terms dominate and almost balance each other. 
The other non-zero volume averages, the production/work terms, are smaller. 

The production of the fluctuating flow is assumed generally to be due to thermal effects 
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E 

D* 

P T 

e 

pF 


8848 

3264 

1639152 

1630896 

0 

filter 

width 

E f 

d* f 

P? 


pF 

1/16 

8712 

3216 

1605264 

1559800 

-37386 

1/8 

8351 

3088 

1520776 

1403512 

-109760 

1/4 

7326 

2704 

1302627 

1080189 

-216349 

1/2 

5261 

1936 

916234 

657181 

-254166 

3/4 

3629 

1280 

632675 

418080 

-209019 


Ec 

D *c 

Pt 

e c 

P c 

1/16 

135 

48 

32961 

67605 

34498 

1/8 

466 

160 

107777 

193472 

85342 

1/4 

1254 

480 

260544 

364362 

102924 

1/2 

2333 

912 

424360 

429592 

2853 

3/4 

2819 

1072 

488818 

406334 

-86228 


Es 

D* s 

P s 

r T 

e 5 

P s 

1/16 

3 

0 

929 

3494 

2540 

1/8 

33 

16 

10600 

33913 

23207 

1/4 

270 

80 

75983 

186347 

110073 

1/2 

1255 

448 

298560 

544125 

244862 

3/4 

2401 

896 

517659 

806482 

287344 


Table 1: Volume Averaged Kinetic Energy Equation Terms 


filter width 

1/16 0.015 

1/8 0.056 

1/4 0.172 

1/2 0.405 

3/4 0.590 


Table 2: Fraction of kinetic energy below the filter width. 
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for the large-scale motion and due to the velocity cascading mechanism for the generation 
of intermediate and small-scale motion. The results of the current simulation are able to 
quantify this notion at least for one Ra. The thermal production occurs at scales larger than 
1/8 based on the observation that Py + P $ , the thermal production below the filter width, 
is reduced to 7% of the total dissipation at this filter width. At the largest filter width, the 
equivalent value is 60%. Thus the thermal production, while concentrated at the largest 
scales, is spread over a significant range. 

The velocity production or cascading energy is shown in Figure 3c. As the filter width 
increases, the velocity production terms in each of the 3 kinetic energy equations must 
approach zero. Assuming that these terms have reached their maximum value near the 
largest filter used in this study, only about 20% of the total dissipation (based on P s at 
A/ = 3/4) is involved in the cascading mechanism. 

In Figures 4 and 5 simple schematics of the energy flows are presented. The conversion 
of energy from one type to another occurs over a range of flow scales and is not generally 
local to a specific spacial region. (The effects of a core and boundary-layer region are not 
distinguished in this volume average analysis.) The average transfer of energy from the 
filtered to the cross to the fluctuation equation is done by the velocity production/ work 
terms. For the flow simulated in this study, energy on the average is extracted from the 
filtered equation and added to the fluctuation equation. The cross velocity production was 
either positive or negative depending on the filter width. 

5 Horizontal Averages 

In this section the variation with x 3 of the horizontal averages of the terms of the set of 
filtered kinetic energy equations (Equations 20, 28 and 36) are examined. Plots of the terms 
in the total kinetic energy equation (Equation 15 and Figure 6a) and the total potential 
energy equation (Equation 49 and Figure 6b) are included for reference. 

The first observation is the existence of boundary-layer regions near the wall where the 
various terms have significantly different behavior from the center layer. In a previous 
analysis of this simulation data, a boundary- layer region of thickness 0.23 next to each wall 
was determined [8]. This agreed with experimental data, in particular that of Carroll[2]. All 
of the terms studied in this analysis exhibited a region, significantly different in character 
from the core region, near the walls with approximately this thickness. (Figures 7 to 9) 
In the core region most of the kinetic energy equation terms are relatively constant with 
the major exceptions being the production/work terms. The production/work terms are 
generally smaller in magnitude than the other terms in the same equation with the same filter 
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width. Their variations with x 3 could be due to the low frequency, patterned flow mentioned 
previously that persisted throughout the simulation. In the boundary layer almost all the 
terms are significant and exhibit large variations in X3. The production terms approach 
zero as they near the walls while the dissipation and diffusion terms (due to the viscous 
component, see below) approach a maximum value at the wall. 

The total diffusion, D, in each of the three equations has a similar variation with X3 
(Figure 7). The diffusion is transferring energy from the core to the boundary layer. The 
diffusion in the filtered equation increases as A^ is decreased. However, as Ay is decreased, 
the diffusion in the fluctuation equation becomes smaller until it is negligible. The diffusion 
in the cross equation is larger than that in the fluctuating equation but also decreases with 
decreasing A f. The diffusion is composed of the stress diffusion plus the two component of 
D * — the pressure component (pressure plus kinetic energy in the filtered equation) and the 
viscous component. Both contribute mainly near the wall, although the pressure component 
has a signifcant core contribution in some cases. Figures 8a and 8b are examples of the 
typical variation with x 3 of these terms, which are similar to the same terms in the total 
kinetic energy equation (Figure 6a). The stress diffusion (Figure 8c) has a shape similar to 
the pressure component and is generally the smallest of the diffusion components. 

An example of the dissipation terms is shown in Figure 9a. There was little difference 
in the shape of the dissipation versus x 3 for the different equations and filter widths. The 
thermal production terms also exhibited this consistency of shape. An example is shown in 
Figure 9b. The magnitudes of these terms for different Aj can be inferred from Figures 3b 
and 3a. One other quantity of interest is the kinetic energy itself. The variation with X3 of 
the three kinetic energies was similar although the fluctuation kinetic energy did not exhibit 
the peak in magnitude in the boundary layer as did the other two (Figures 9c and 3d). 
All three showed no significant variation in shape with A f. The variation with time of the 
kinetic energies can be inferred from Figure lb. 

6 Kinetic Energy Equations 

All the terms of the filtered kinetic energy equation are plotted together in Figure 10 for 
3 filter widths. The dissipation, e, is positive when it extracts energy, which is opposite to 
the other terms — due to writing — e for the dissipation in the kinetic energy equations. The 
terms have similar variations in X3 to those of the terms in the total kinetic energy equation 
(Figure 6a), except that the velocity production terms are non-zero although they are small. 
The volume averages for these same terms except for the diffusion which is zero are plotted 
versus A/ in Figures 11a and lib. While the thermal production balances the dissipation 
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globally (volume averaged), this is not true locally (horizontally averaged), particularly in 
the core. For small A / the dissipation in the core increases but still does not balance the 
production. 

For a large A/, the filtered equation describes the energy flow of the large scales. From 
Figures 12a and 12b the energy extracted from the large-scale flow by the production terms 
is mainly in the boundary layer. For the smaller A / in Figure 12c the fluctuations due to the 
scales smaller than the filter width are extracting energy more uniformly across the entire 
fluid layer. Another observation is the change in importance of the three production terms 
shown in Figure lib. At small A } , and P[ ! are dominant. At the larger A/ where Cij 

and ~Lij approach zero, P R becomes largest — although this term also must approach zero 
for very large A/. 

The cross kinetic energy equation, Equation 28, contains 9 pro Auction/ work terms. For 
horizontal and volume averages, a simplified form reduces the number of terms to 3 (Equa- 
tion 59). Since no particular physical significance can be associated with either form, the 
data will be plotted for the reduced form for simplicity. 

The cross kinetic energy equation terms are shown in Figures 13 and 14. For larger A/, 
the variation of the terms with x 3 is similar to that in the filtered equation (Figure 10). 
One difference is that the dissipation is more significant in the core for the cross equation 
than for the filtered equation. For smaller A/, the velocity production is almost equal to the 
thermal production; whereas, in the filtered equation the thermal production dominates the 
velocity production for all A/. This is due to the large magnitude of the velocity production 
component, P%f • 

Similar plots for the terms in the fluctuation kinetic energy equation are shown in Fig- 
ures 15 to 17. To simplify the plots, the two velocity production terms involving C\j have 
been added together. The same was done with the R^j terms. In each case, one of the two 
terms is significantly smaller. 

Pc » p c ( 60 ) 

Pr « Pr ( 61 ) 

The data for the fluctuation equation is similar to that for the cross and filtered equation. 
However, the boundary layer and core regions are not as well defined as in the other two 
equations. On close inspection there appears to be a region in the core merging with the 
boundary layer — 0.2 <C | X3 1 <C 0.4 — where the velocity production and the diffusion have 
increased importance. Another difference from the other two equations is that the thermal 
production does not approximately balance dissipation even on a volume averaged basis. 
The velocity production is even slightly larger than the thermal production for small A/. 
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The magnitude of P s is not mainly due to one component, as for the cross equation, but 
the “L, C and R” components vary in importance as a function of A / (Figure 17b). 

7 Subgrid Scale Modeling 


Calculation of the terms in the kinetic energy equation using the direct simulation data can 
be used to study subgrid models for LES. For example, the production terms for Hij+Uij can 
be calculated directly (refered to as “exact”) and compared to the same term calculated with 
a model for + Cij. The “exact” and modeled production terms can then be compared 
using a correlation analysis[l, 4], However, a direct comparison of some average of the 
production terms (such as an average over homogeneous directions) can be used also to 
evaluate the models. Since the selection of models and model constants is often based on the 
correct prediction of energy transfer between different flow scales in many LES studies [1], the 
correct prediction of the magnitudes of the production terms should be an important model 
criterion. In particular, Pr+c = Pr + Pc * > which extracts energy from the filtered flow 
field (Equation 20), should provide the most important criterion. The analysis in this study 
is only the first step in evaluating subgrid models. Conclusions drawn from this analysis of 
DNS data which suggest changes in subgrid modeling must be verified in an actual LES. 

The main thrust of this study was to evaluate the behavior of commonly used subgrid 
models for natural convection flows. To the authors knowledge, most studies of subgrid 
modeling using DNS data have been done using isothermal flows. The model comparisons 
were made using a range of filter widths from 0.0625 to 0.75 . These filter widths were chosen 
so that the split of kinetic energy between the filtered velocity field and the subgrid velocity 
field would include the range found in typical LES studies. The Smagorinsky model and the 
scale-similarity model were the models tested since these are the two most frequently used 
models. The Smagorinsky model is [7]: 


M - = -K~5„ 


(62) 


* = 

\/2 

^ = (SijSij)* 

A s = minimum (A f or distance to nearest wall ) 
Ck = model constant 

The scale-similarity model is [1]: 


( 63 ) 
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filter width 


Ck 

1/16 0.047 

1/8 0.055 

1/4 0.069 

1/2 0.114 

3/4 0.175 


Table 3: Smagorinsky model constant, C K , adjusted so that the model predicts the “exact” 
value of Pr+c eac h filter width. 


u' x u' } - \SijU' k u‘ k = (UiUj - UiUj) - - UnUn) 


(64) 


The value of Ck for the Smagorinsky model is one of the results of the current analysis. 
Generally, Ck — 0.21 has been used in natural convection simulations [7], but this value 
assumes that the grid size, A, is used as the length scale, A,. Studies for homogeneous 
turbulence have shown that Ck should be reduced by a factor of 2 or 3 when A f is used 
instead. Therefore, a value between 0.07 to 0.14 should be assumed as a target value of 
Cr- If the constant of the Smagorinsky model is adjusted to match the same volume- 
averaged production as the “exact” term, the value of Ck varies as shown in Table 3. The 
Smagorinsky constant can also be allowed to vary with x 3 . If Ck is adjusted at each 13 to 
give the simulation results, the modified Ck varies with x 3 as shown in Figure 18. A value of 
0.05 gives best agreement in the center region and the value of Ck varies significantly closer 
to the walls. Ck in the side regions increases with filter width. 

< Pj/ +C > is plotted in Figure 19 for the simulation data (labeled “exact”) and for the 
two models mentioned above. < Pr+c ^ f° r f be Smagorinsky model was calculated at the 
different filter widths with the value of C K in Table 3. Since the scale similarity model has 
no adjustable constant, {Pr+c) calculated using this model is different from the other two 
in the horizontal average plots. {Pr+c} versus filter width is plotted in Figure 20 where the 
same value of C K (0.07) was used for all filter widths. The data in Figure 19 gives a more 
local comparison of the models to the “exact”. At the small filter widths, the scale- similarity 
model gives good agreement with the “exact” value as x 3 (the inhomogenous direction) is 
varied. For larger filter widths, the agreement is not as good and the scale similarity model 
even predicts negative values for < Pr +c > in the core region while the “exact” value 
remains positive. The Smagorinsky model shows no agreement with the variation of spacial 
distance in an inhomogenous direction and can only be providing a global agreement in 
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filter width Smagorinsky scale-similarity 


1/16 

1/8 

1/4 

1/2 

3/4 


0.07 

0.12 

0.19 

0.19 

0.14 


0.99 

0.97 

0.89 

0.74 

0.68 


Table 4: Volume average of the correlation of Pr+ c calculated directly with Pr+q calculated 
using a model for Rij and Uij. 


energy dissipation as has generally been believed. For comparison of these results with other 
calculations, the fraction of energy below the filter width shown in Table 2 may be a better 
parameter than the value of the filter width itself. 

Previous studies have estimated that the scale- similarity model has a very small contri- 
bution to the subgrid turbulent production and dissipation [1]. The current results show 
that the fraction of the “exact” production energy predicted by the scale-similarity model 
varies significantly with A f (Figure 20). For the smallest A j tested (too small to be of use 
in a LES), the scale-similarity model prediction of < Pr+c > i s almost equal to the “exact” 
value (Figure 19). As A / increases, this fraction predicted decreases to 12% for A / = 0.75. 
As mentioned above, the value of Cr that best matches {Pr+c} °f the Smagorinsky model 
to the “exact” value also varies with A f. In most LES studies the selection of a model and 
the model constants is reported without reference to the fraction of energy in the subgrid 
field. Based on these results it appears that the choice of models and model constants should 
be a function of the energy fraction in the subgrid field. 

The correlations of the two models with the “exact” term based on the simulation data 
for < Pj/ +C > are shown in Figure 21. The horizontal correlations are formed as follows: 


corr(x 3 ) = 


< (a- < a >){b- < b>)> 

,J< (a- < a >) 2 >\J< (< b— < b >) a > 


(65) 


where a and b are functions of X 2 and X3. The correlations of Pr+c an< ^ ( no ^ 

shown) have roughly the same magnitudes. The correlation of Wr+ c is the same “scalar 
level” comparison in the study done by Bardina et. al.[l], except that the correlations were 
done in horizontal planes in the present study. The vertical average of the results in Figure 21 
are shown in Table 4. The correlation levels in the current study show a significantly higher 
correlation for the scale-similarity model than for the Smagorinsky model, similar to that in 
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[1]. In [1] the correlation coefficients were 0.50 and 0.58 for the scale similarity model and 
0.36 and 0.05 for the Smagorinsky model. The flows used in that study were homogenous 
isotropic turbulence and homogenous shear turbulence. The current correlation results are 
in general agreement with [1]. The scale similarity model gives better local agreement with 
the “exact” results than the Smagorinsky model. 

8 Concluding Remarks 

A formal study of the set of kinetic energy equations derived from filtering the original 
momentum equations has provided insight into the physics of the current flow. In particular, 
the variation of the thermal and velocity production terms with scale (or A/) has been 
quantified. When one deals with the large quantity of detailed data that exist in turbulent 
flows some averaging of the details is necessary to study the flow. Fourier spectral analysis 
provides information on the spacially averaged energy flow. The filter analysis presented here 
provides similar information but with a spacially local average. Also, the current study has 
shown that the diffusion terms are not negligible, especially near the walls. The Smagorinsky 
models assumes production equals dissipation, neglecting diffusion. 

The extension of the calculation of the terms of the filtered set of kinetic energy equations 
to evaluate subgrid models is natural. One difference in this study from previous “a priori” 
studied is that the fine scale calculations of the direct simulation were not projected onto a 
coarse grid before calculating the model results. This was done to decouple the effect of the 
subgrid filtering from the filtering that results from projecting onto a coarse grid. 

The validity of the “a priori” studies to investigate LES models using DNS data is an 
open issue. Possibly the approach taken in this study when applied to a wider range of flows 
and then compared in actual LES studies can help clarify the process of determining subgrid 
models. In particular, a better understanding of the effect of combining the Smagorinsky 
model with the scale similarity model (called the linear combination model) is needed. In 
the study where the linear combination model was proposed [1], it was suggested that the 
scale similarity part contributed little (about 5%) to the energy dissipation. The current 
results suggest that the amount varies with A / and is generally a larger percentage. 

9 Appendix A - The Rayleigh- Benard problem 

The Rayleigh- Benard problem is a simple geometry, laboratory-type problem used to study 
natural convection. The problem entails the study of the fluid motion and thermal convection 
of a rectangular fluid layer which is typically heated from below. The layer is typically 
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thin, has no mean pressure gradient in the horizontal directions, and has uniform boundary 
conditions in the horizontal directions. 

Wall boundary conditions are assumed at the upper and lower fluid edges. Periodic 
boundary conditions are assumed in the horizontal directions. The vertical coordinate, x 3 , 
was defined as zero halfway between the two walls for this study. 

1. Variables: 

• fluid properties (all constant) - 

v kinematic viscosity 

a thermal diffusivity 

/ 3 coefficient of thermal expansion 

p 0 reference fluid density 
g acceleration of gravity 

• geometry constants - 

h distance between the two horizontal walls 

• dependent variables - 

Ui fluid velocity, i = 1,2,3 
T„ fluid temperature 
P a fluid pressure 
p fluid density 

• independent variables - 

X{ spacial coordinates, i = 1,2,3 
t time 

• flow constants - 

T a temperature at lower wall 

AT temperature difference between upper and lower wall 

2. Boundary conditions: 

Uj = 0 x 3 = ±h/2 

T a - T 0 x 3 = -h/2 

T 0 = T 0 — AT x 3 = h/2 

AT >0 

3. Boussinesq assumption: 


P9 = P<>9 - p 0 gP(T a - T 0 ) 


4. Removal of static temperature and pressure gradient 


( 66 ) 


T 0 (x,f) = T r (x,t) + T 0 


(67) 
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( 68 ) 


AT 

T r {x,t) = T(x,t) - ~2j^( x 3 + h ) 


d?a 

dxi 

dP a 


dP 

dxi 

dP 


+ Xf, i — 1,2 

AT, 


+ 


Po9 (l + P- 2 h( x i + M) 


dx 3 dx 3 
5. Non-dimensionalize by a, AT, p a and h 


Ui = 


P = 


T = 


a/h 

P 

p 0 a 2 /h 2 

T 

AT 




va 


PrEE- 

a 


Xi = 


t = 


h 2 /a 

6. Governing Equations 
dui 

TF 1 ^ 0 

OXi 


dui A dxii dP _ dSij _ _ __ 

"57 + + Pr ~x^ 1 + PrRaT6 i3 

at oxj OXi oxj 

d f * Qf Q2 f * 

H + Ui W i = dE? + V3 

q _ OUj 

ij ~ dxj + dx, 


(69) 

(70) 

(71) 

(72) 

(73) 

(74) 

(75) 

(76) 

(77) 


(78) 

(79) 

(80) 
(81) 


In the main body of this paper the non-dimensional variables and equations will be used. 
The * notation will be removed for convenience. 
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( 82 ) 


10 Appendix B - Production Terms 


pft — ./nr. . 

r c — 

H’ a (83) 

# = (84) 

Pc' = - (85) 

P/S' = -jS/Jla (86) 

P" = -iSaBs (87) 

P% = (88) 

PH = —SlA ( 89 ) 

Pl< = -±S. y I y (90) 

P£' = (91) 

= -isitfiyiiy) (92) 

P&' = -1^(535) (93) 
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Figure 3: Volume- averaged terms in the set of filtered energy equations 
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Figure 4: Volume- averaged energy flow diagram for filter width 0.25 . All 
data has been divided by 10®. 
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(a) kinetic energy equation 
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(b) potential energy equation 

Figure 6: Horizontally-averaged terms of the total kinetic and potential en- 
ergy equation 


3 





32 







Figure 8: Horizontally-averaged components of the diffusion term for A / 
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Figure 9: Horizontally-averaged terms in the set of filtered kinetic energy 
equations for A / = 0.5 
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Figure 10: Horizontally-averaged terms in the filtered kinetic energy equation 
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Figure 11: Volume-averaged terms in the filtered kinetic energy equation 
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Figure 12: Horizontally-averaged components of the velocity production term 
in the filtered kinetic energy equation 
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Figure 13: Horizontally-averaged terms in the cross kinetic energy equation 
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Figure 15: Horizontally-averaged terms in the fluctuation kinetic energy 
equation 
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Figure 16: Horizontally-averaged components of the velocity production term 
in the fluctuation kinetic energy equation 
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Figure 17: Volume-averaged terms in the fluctuation kinetic energy equation 
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igure 18: Ck versus x 3 for modified Smagorinsky model 
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Figure 19: Horizontal average of Pj/ +C - simulation and 2 models 
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Figure 21: Correlation of 2 models with simulation results for PtS, 
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